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^p Abstract. We introduce a method for calculating the stationary state of a translation 

invariant array of weakly coupled cavities in the presence of dissipation and coherent 
OQ as well as incoherent drives. Instead of computing the full density matrix our method 

directly calculates the correlation functions which are relevant for obtaining all local 
j— | quantities of interest. It considers an expansion of the correlation functions and 

f^h their equations of motion in powers of the photon tunneling rate between adjacent 

cavities, leading to an exact second order solution for any number of cavities. Our 

method provides a controllable approximation for weak tunneling rates applicable to 

C3 the strongly correlated regime that is dominated by nonlinearities in the cavities and 

thus of high interest. 
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1. Introduction 

The study of light matter interactions that are enhanced by confining light fields in 
electromagnetic cavities has been a thriving discipline of Quantum Optics throughout 
the last decades. Particularly with the advent of realizations of the so called strong 
coupling regime where the strength of the interaction between a photon and a 
quantum emitter exceeds the decay mechanisms for photons and emitters, experimental 
investigations of the coherent interactions between single emitters and individual 
photons became possible pp. 

In recent years, a new direction of cavity quantum electrodynamics (cavity-QED) 
has developed, in which multiple cavities that are coupled via the exchange of photons 
are considered. Such setups are particularly intriguing if the cavities are connected 
forming an array and the strong coupling regime is achieved in each cavity of the 
array. These devices would then give rise to quantum many-body systems of strongly 
interacting photons and polaritons [HI HI [3]. As an alternative to a cavity array, one 
may also consider optical fibers that couple to nearby atoms [6j [7j or even clouds of 
Rydberg atoms that are optically thick in free space [SJ. Both these systems avoid 
the need to build mutually resonant cavities, which is possible [9] but can be rather 
challenging in the optical range. For microwave photons it is however perfectly feasible 
to build large arrays of mutually resonant cavities on one chip in an architecture known 
as circuit-QED HDUTT]. 

For strongly interacting polaritons and photons in coupled arrays of micro-cavities 
and optical fibers, possibilities to observe equilibrium phenomena, such as a Mott 
insulator [HJIHE] or a Tonks-Girardeau gas [HlEj, have mostly been addressed so far and 
the development has been summarized in the reviews [12l [13J [10] . In every experiment 
that involves light-matter interactions, some photons will however inevitably be lost 
from the structure due to imperfect light confinement or emitter relaxation so that 
thermal states are no longer an appropriate description of the system. To compensate 
for such losses, coupled cavity arrays are thus most naturally studied in a regime where a 
coherent or incoherent input continuously replaces the dissipated excitations. This mode 
of operation eventually gives rise to a driven dissipative regime, where the dynamical 
balance of loading and loss processes leads to the emergence of stationary states. Yet 
the properties of stationary states of driven dissipative systems are only explored to 
a much lesser degree than the properties of thermal equilibrium states. For coupled 
cavities, small arrays have been considered exactly [H] and mean field approaches for 
larger arrays have been employed [TS]. Moreover numerical studies found signatures 
for crystallization [16] and photon solids were predicted for arrays with cross Kerr 
interactions [T7] . 

As driven dissipative quantum many-body systems have to date only barely been 
explored, there is a need for technical tools for their efficient description. Here we 
introduce a perturbative technique for the calculation of the physical properties of 
stationary states of driven dissipative cavity arrays. Our approach assumes a large, 
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translation invariant cavity array where the quantum states of all cavities are identical. 
Instead of computing the full density matrix, it directly calculates the correlation 
functions which are relevant for obtaining all local quantities of interest. A power 
expansion of the formal exact solution to second order in the photon tunneling rate 
between adjacent cavities, provides semi-analytical results for any number of cavities. 

2. Exact solution in the steady state 

We consider an array of cavities that are coupled via mutual photon tunneling and 
where each cavity is doped with a Kerr nonlinear medium that generates a strong 
photon-photon interactions. In a frame that rotates at the frequency of the coherent 
drive lasers, this system is described by a Bose-Hubbard Hamiltonian with local coherent 
drives (h — 1), 

U 



H N = ^2 A Q l a J + -w a l a l a i a i + ^( a l + a i 



=i 

iV-l 

E 

=i 



+ J y2(a\a i+1 + a\ +1 ai) + J((J N a-y + a\a N ) (1) 



where we assumed periodic boundary conditions so that N cavities form a circle coupling 
to nearest left and right neighbors. Here, A = u a — wl is the detuning of the cavity 
resonance frequency u a from the laser frequency u^, U is the strength of the Kerr 
nonlinearity, J the rate of photon tunneling between the cavities and Q the drive 
amplitude of a coherent laser drive. We consider the local energy scales orders of 
magnitude smaller than the cavity frequency, that is, U, J, A, Q <C uj a , so that the 
rotating wave approximation can be applied. Adding decay and incoherent pumping of 
the modes, the total Liouvillian that describes the dynamics of this system reads 

N 

d t p = i[p, H N ] + - y^(2o i /pq| - a\aip - /Sajoj) 
i=i 

P N 
+ ^ ^(2o}pOi - a t a\p - pojoj) , (2) 

i=i 

where p denotes the total density matrix of the cavity array. 7 is the rate of photon 

decay and P the rate at which photons are incoherently pumped into the device. We 

are interested in the steady state under a continuous excitation of each of the cavities, 

represented by the reduced density matrix p of a single cavity. As all cavities have 

exactly the same properties and dynamics, we find pi — P2 — ■ ■ ■ — Pn — P, where pj is 

the reduced density matrix of cavity number j. 

Here, instead of obtaining the full density matrix of the array in the steady state 

{d t p = 0) and then tracing out all cavities but one, to obtain p, we compute the steady 

state properties directly in the form of the local mean values of all possible operators 

defining the system. These can be written as 

(4 m a?) = Tr(parX) (3) 
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where m and n are integers. Let us call O the set of operators the averages of 
which correspond to the correlators required to describe the full iV-cavity system, 
i.e., O includes all the sought observables as well as operators which couple to them 
through the equations of motion, (a\ m a"a 2 >1 a v 2 . . . a^a N ). To simplify notation, we 
express this general correlator as {{m, n}, {/i, z/} ... {a, (3}}. Naturally, in the case 
of an anharmonic mode or in the presence of any anharmonicity, one must choose 
an appropriate truncation in the number of excitations, i.e. n,m < n max , that also 
truncates the number of correlations in this set. For a given driving intensity, the 
truncation must be high enough to yield converged, accurate results. 

From the master equation ^ one can obtain the set of coupled equations for the 
full set of operators O, which we write in the matrix form, 

d t v = Mv + I } (4) 

with all correlators forming the vector v, i.e. v = ((ai), (a[), . . . , (aiG^), {a\a^) ■ ■ .), 
where the exponent T denotes the transpose. The coefficient matrix M and vector / 
are derived from the master equation in a systematic way [181 HH] , as we explicitly show 
in the Appendix. The solution of Eq. (ji!) is completely equivalent to computing the 
correlators as in Eq. pL from the density matrix obtained by solving Eq. (pi). The 
exact solution in the steady state (if it is unique and exists) simply reads 

v = -M- 1 ! (5) 

We can considerably reduce the number of operators to a minimal set by making 
use of the translational symmetry in the ID chain (circle). That is, all correlators that 
are left or right circular rotations of the elements in {{m, n}, {//, v\ . . . {a, (3}}, such as 
{{a, (3}, {m, n}, {/i, u} . . .} or {{a, (3}, . . . , {/i, u}, {m, n}}, are redundant because they 
are exactly the same. There is, therefore, a maximum of 2N representations of the same 
correlator (less if some of the pairs are {0,0} or mutually equal). We can choose an 
arbitrary rule to systematically keep only one representative of such set of redundant 
correlators, for instance, we choose the ones where: 

(i) The nonzero sets are always to the left and as cluttered together as possible, such 
as {{4, 2}, {3,0}, {0,0}...} 

(ii) The largest sum of indexes is most to the left: m + n > \x-\-v > ... > a + (3, such as 

{{3, 3}, {1, 3}, {0, 2}, {1, 0}}. Together with the previous rule, this may give things 

like {{1,3}, {3, 3}, {1,0}, {0,0},...}. 
(iii) If there are two pairs with an equal sum, the one with the largest first index 

is left-most m > fx. Together with the previous rule, this may give things like 

{{4,2},{3,3},{0,0},...}. 

With this, the vectorial Hilbert space is significantly reduced, for instance from 6559 to 
1033 for N = A and n max = 2. From v we thus extract a new vector v that only contains 
the minimal set of correlators. The equation of motion for v is then written in terms 
of a new matrix, M , reconstructed by removing the redundant rows and summing the 
coefficients of the redundant columns of M, and a new vector, I removing the redundant 
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rows of /. The exact stationary state solution of the reduced system of equations is now 
given by 

V = -M~ 1 I. (6) 

2.1. One cavity, N — 1, and the uncoupled limit 

In the uncoupled limit, J = 0, the calculation can be reduced to a single cavity. The 
operators for each cavity, {a™ a™), are the same for all i — 1, . . . , N, so let us denote 
them by a common name, (a) m a n ) = {{m,n}}, without any index. We can define the 
vector v a of all possible individual cavity operators, 

Wa T =((a),(at),(ata),...) . (7) 

In the case of one cavity or many which are uncoupled, the full ensemble vector v reduces 
trivially to v = v a . The equation of motion for this system reads d t v a = M a v a + I a and 
its stationary solution is, 

v a = -M~ 1 I a . (8) 

2.2. N < 4 cavities 

Four is the minimum number of identical cavities needed to obtain a general solution 
to second order in J that is valid for any N. The reason is that four coupled systems 
imply a qualitative change, as compared to two or three, as each of them is no longer 
in contact with all the others. This represents the general case to second order in J. 

In this case, the vector of correlators v not only contains the subset v a with 
{{m,n}, {0, 0}, {0, 0}, {0,0}}, but also another four new subsets that include cross 
correlations with two, three or four cavities. The first one, which we call i>&, includes 
correlators with two cavities, (a i 0^02 a%) = (oj ajoj oj) = . . . = (a^a^b^b 11 ) = 
{{m, n}, {fi, u}, {0, 0}, {0, 0}}. We denote with b the photon annihilation operator for 
the second cavity, and apply the rules to extract the minimal set of operators (m + n > 
\x + v, with m > fj, in case of degeneracy of the sum). For example, for a truncation in 
each cavity system with n max = 2 photons, we thus have a dimension of 8 for v a and 36 
for vj,. A third subset, v c , includes correlators with three cavities, (aj" 1 OiOj a^a^ a|) = 
(4 m a£ at /X a 3 Pa 3> = • • • = (a^ m a n b^b u c^c q ) = {{m, n}, {pi, u}, {p, q}, {0, 0}}. Similarly 
to b, we denote with c the photon annihilation operator in the third cavity, where 
we have applied the circular rules to obtain the minimal set of operators. Finally, 
specifically to iV = 4, we need to consider Vd, which includes correlators with 4 cavities 
(a^a^a^a^a^alal'al) = ... = (a^ m a n b^b u c^c q d^ s d t ) = {{m,n}, {//, u}, {p, q}, {t, s}}, 
d is the photon annihilation operator in the fourth cavity, and v e which includes 
operators of two cavities at alternate positions, (a| m a"a3 a|) = (a 2 m a^a/ a|) = . . . = 
{a) m a n c^c q ) = {{m, n}, {0, 0}, {p, q}, {0, 0}}. 

We can rewrite Eq. (Ek in terms of these subsets of correlators, each with a different 
dimension, as five coupled matrix equations, 

d t v a = {M a + US a )v a + I a + iJR ab v b , (9) 
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d t v b = (Mb + iJS b )v b + (B ba + iJR ba )v a + iJR bc v c , (10) 

d t v c = (M c + iJS c )v c + (B cb + iJR cb )v b + iJR cd v d + (B ce + i.JR ce )v e ,(11) 
d t v d = (M d + iJS d )v d + (B dc + i JR dc )v c , (12) 

d t v e = (M e + US e )v e + (B ea + UR ea )v a + UR ec v c . (13) 

Here, we have separated the effect of the hopping J into the self-renormalization matrices 
S, and the linking matrices R, that only contain integer numbers. The vector I a and 
matrices B, contain only the driving parameters Q and P (coherent or incoherent). 
Other internal parameters such as A, 7 and U enter in the matrices M. 

One can solve these equations recurrently in the steady state, from bottom to top. 
This may be useful for a small number of cavities where the expressions are simple. For 
instance, for N = 2, the system reduces to v a and v b with Eqs. (J9])-(10), and solutions: 

(14) 
(15) 



V a = - (M a + iJS a + iJRahFha)' 1 I a , 

v b = F ba v a , 

where F ba = — (M b + iJS b )~ (B ba + iJR ba )v a . Similarly, for iV = 3 we have: 

v a = - (M a + i JS a + iJR ab F ba y l I a , (16) 

v b = F ba v a , and v c = F cb v b , (17) 

where F ba = - (M b + iJS b + iJR bc F cb y l (B ba +iJR ba ) and F cb = - (M c + iJS^ 1 (B cb + 
iJR cb ). This recursive procedure is possible in principle for any N, although, in gen- 
eral not very practical, given that the exact solution can also be obtained by simply 
inverting one matrix, M, as in Eq. ([6]). Anyhow, obtaining the exact solution becomes 
exceedingly cumbersome for a large number of cavities, iV 3> 1. In the following we 
therefore concentrate on finding an approximate solution. 



3. Approximated solution to second order in J 

In order to find an approximate semi-analytical expression for the steady state of a 



cavity, v a , we expand both the correlators in v and the set of equations <\9\)-{ 13 ) in powers 



of J up to second order. More precisely, second order for v a , requires for consistency 
the following lower orders in the other subsets: 



v a 
v b 

Vc 

I'd 

V e 



vi 0) + Jv a l) + JV + 
vi 0) + Jv^ + . . . , 

(0) 



V 



(I 



+ ... 



v^ + 



(18) 
(19) 
(20) 
(21) 
(22) 



The expanded equations read in these terms: 



d t v a = 



M a v a 0) + I a ]+J \M a v a l) + iS a v a 0) + iR a , 



bV b 



(0) 



+ J 1 



M a v a 2) + iSaV^ + iR ab v[ 



(1) 



+ 



(23) 
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d t v b = 



J 



M b vf ] + B ba v^ 
M b vl 1] +iS b vl 0) 



,(o) 



B ba v^ + iRtovW + iR bc vf +.... (24) 



(o) 



We truncate the equations at this point since the solutions obtained by setting each 
square bracket to zero, 



>i° } = - M~H a , „<°> = -M^B ba v^ 



(1)_ 

a 


- K" 1 


'i-S^f+zi^f 


5 


f = 


-M," 1 


'i^^ + S^+i^f 


(2) = 

a 


-M- 1 


iSavP+iRabvP 


5 



zi?^ ) 



(25) 
(26) 

(27) 

(28) 



,(o) 



(0) 



ultimately depend on v c only. Obtaining v c from the equations would in turn require 
the knowledge of v e but this is not needed given that the zero order is simply the 
uncoupled limit, that is, products of the solutions for N = 1 as in Eq. (J8J). For instance, 
the uncorrelated solution for v h , corresponding to (a) m a n b^b u )^ = (a^ m a n )(a^ fl a u ), can 



V a s\. b V a , 



where X b is 



be directly obtained through the product of twice v i , as v j, 
the corresponding mixing matrix obtained by inspection. This is completely equivalent 
to the linear algebra solution of Eq. (l25). The same applies for v c but with two mixing 
matrices, vf' = i4° X cl Va X c2 Va ■ 

These solutions are valid for N > 4, since adding more cavities to the circle does 
not produce any structural qualitative change to second order in J. The approximation 
is better the larger the N. It is formally the same for N = 2 and 3 (setting Vc — for 
N = 2), but differs quantitatively to first and second order, respectively due to different 
coefficients in the equations. 



4. Comparison between the exact and approximated results 

Since the approximated solutions are single valued, this method cannot reproduce 
regimes where several steady states are compatible for the individual system 
(corresponding to different steady states of the ensemble) or any other instability regions 
like lasing. Its perturbative nature allows it only to describe regimes where the coupling 
is smaller than the effective decoherence or driving. More precisely, the weak coupling 
regime, where new collective eigen-modes are not required to describe the ensemble 
dynamics. 

We illustrate the interest of this method by comparing in Fig. [l] the exact solution 
for iV = 4 (in solid black) with the approximated ones, valid for N > 4. We have chosen 
two quantities of interest, the mean cavity population, n a = (a^a), and its second order 
coherence function at zero delay, g^ 2 \0) = (a) a) aa) / 'n\. Figs, hi a) and (b) show that 
both are well approximated by the second order solution (in dashed red) as long as 
J < n, 7. Lower order approximations (in blue and green) deviate from the exact 
solution at even lower J. In Figs. [11(c) and (d), we fix J = O.37 < fi = O.57, and scan 
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Figure 1. Mean cavity population, n a , and second order coherence function, </ 2 )(0). 
The exact solution for N — 4 cavities (solid black) is compared to second (dashed red) , 
first (dashed blue) and zero (dotted green) order approximations, valid for N > 4. In 
(a) and (b) we fix fl = 0.77 and A = and vary the photon tunneling rate J. In (c) 
and (d) we fix fi = O.57 and J = O.37 and vary the laser frequency wl- In the inset 
of (d), we have magnified the region around 0. Other parameters: U = 67, P = 0, 



the system resonances by tuning the laser frequency. In this case, we observe that the 
second order approximation remains very close to the exact solution for all frequencies 
while the first and zeroth order deviate from it, notably, close to the various cavity 
resonances at 0, U/2, U (marked with vertical lines). The first order approximation 
breaks down near the one-photon resonance at as evidenced by the negative value of 
g( 2 )(0) in the zoom-inset of Fig. Hid). 

Due to the form of the dissipation terms in Eq. (J2J) we cannot illustrate the method 
in the absence of a coherent drive, only under the action of an incoherent pumping. An 
incoherent pump bringing the system into a steady state is equivalent to letting each 
cavity interact with a thermal bath, where P = 71t7o and 7 = (1 + ^t)7o- Here 70 is the 
decay rate at zero temperature and nx the occupation number of the bath at temperature 
T, which are identical for all cavities in our considerations. Hence for the case of purely 
dissipative dynamics with Hjy = in Eq. d2l), the steady state p^ is a product of 
thermal states at temperature T for each cavity, that is, p+j. = Z~ 1 ^2 n e~ ulan ^ kBT \ 
where ks is Boltzmann constant, n the total number of excitations in the system and 
Z the partition sum. For Vt = the Hamiltonian in Eq. ([!]) conserves the number 
of excitations in the system and thus [H N ,p^] = so that p^ is the steady state 
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Figure 2. (a) Mean cavity population, n a , and (b) second order coherence function, 
</ 2 )(0), as a function of the thermal bath occupation. The exact solution for N = 4 
cavities (solid black) is compared to second (dashed red), first (dashed blue) and 
zero (dotted green) order approximations, valid for N > 4. Parameters are chosen 
to maximize the cavity population in Fig. [He): U = 670, Q = O.570, J = O.370, 
wl = w a + 0.257o and n max = 2. 



even in the presence of unitary dynamics generated by H^, independently of the value 
of J. Of course the rotating wave approximation that has been applied to derive the 
Hamiltonian ([!]) is only valid for (/, J C w a . In this regime a thermal bath as described 
by the dissipation terms in Eq. (pi) is "blind" to the energy scales U, J and all cavities 
will eventually be in thermal equilibrium with their bath, n a = n^ and g^ 2 \0) = 2, 
regardless of the hopping J and therefore other neighboring cavities. The dynamics and 
spectrum of emission (out of the scope of the present study) do, however, depend on the 
microscopic properties of the cavities. For example, larger J and 70 would accelerate 
the thermalization of the cavity array. 

If, on the other hand, a coherent and an incoherent drive are both present, the 
steady state becomes nontrivial. Moreover, our approach is suited for exploring this 
experimentally relevant scenario that describes coherently driven cavities in the presence 
of thermal background radiation. Fig. L^shows n a and g^ 2 ' (0) as a function of the thermal 
bath occupation number %, which increases with increasing temperature, for the case 
of maximum cavity population in Fig. flTc). Both the exact and the approximated 
solutions converge at high temperature to the thermalized steady state (n a = ut and 
g( 2 \0) =2). At low temperatures, where the emission is antibunched and each cavity 
behaves like a single-photon emitter, first and zero order approximations differ strongly 
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from the exact solution, especially when computing n a , while, again, the second order 
approximation follows it quite smoothly. 

It is interesting to note that the steady state cavity spectrum of emission could also 
be obtained from this approach without recurring to the quantum regression theorem 
and, therefore, to deriving any time dynamics. The alternative to such complications 
is to look into the steady state occupation, as a function of its natural frequency, of 
a mode, that is weakly coupled to the cavity array and which plays the role of the 
detector. We showed with coworkers the equivalence between this quantity and the 
power spectrum [T9] . 

5. Conclusions 

We have presented a method to solve for the steady state of coupled cavities in a circular 
ID array with translational symmetry, to second order in the photon tunneling rate, J. 
This method can be generalized to any set of identical weakly coupled systems, being 
in a ID, 2D or 3D arrangement. 

We consider any type of driving of the cavities (coherent or incoherent), dissipation 
and nonlinearities. We first derive the equations of motion for a minimal set of relevant 
correlators, v a , and then perform a power expansion of both the equations and the 
solutions to obtain semi-analytical expressions for v a ~ v a + Jva + J 2 v a ■ The 
approximated solution is invariant for N > 4 cavities due to the nearest neighbor nature 
of the coupling. We have finally illustrated the performance of our method with an 
example of four weakly coupled cavities under a coherent drive and temperature. 
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Appendix: Equations for the correlators 

In this appendix we provide the matrix M and vector I appearing in Eq. (kj), which are 
the starting point of the described procedure. Let us consider operators for only two 
adjacent cavities, (ffl}" 1 ajaj flj), as the general case of an array is an straight-forward 
generalization. Then, we can obtain their equations of motion from the full master 
equation as 

d t ^a\afa v 2 ) = Tr ($p afX4X) = V R „,,„,,.,„ (a\ k a\a\ a a^) . (29) 

k, p, a, j3 
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The corresponding elements in R are given by [18] : 

7 U 

R m ,n,v,v = iA(m — n) — — (m + n) + i — [m(m — 1) — n(n — 1)] 

m, n, fi, v Zs Zs 

+ iA(fM -u)-fa + v) + i|[^ - 1) - v{y - 1)] , (30) 

R m,n,n,v = Pmn , -R m,n, M ,„ — P [W , (31) 

m — 1, n — 1, jU, v m, n, /j, — 1, 1/ — 1 

-R m,^,* —iU(m — n), R m , n ,^ v = iU(/j, — v) , (32) 

m -J- 1, n -}- 1) £*, f m,n,/^ + l,i/ + l 

-it m,n,fi,v — %\ LTfl , XL m,n,^,u — Zi Z^X , V^^J 

m — 1, n, /j,, v in.n, /j, — l,v 

-fi/ m,ri,i±,v tli/i ; -fi m,n, jj,,u v<\Lls ; ^Ot:J 

m , n — 1, /j,, u 7n, 7i, /j,, i/ — 1 

-fi m, n, (j,, v to 1 1 1 j -Lb 7n, n, j_i, v LU fJj ^ \O0 J 

m — l,n,n,-\-l,u m + l,n,^i — 1,1/ 

-fi m, n, /j,, v Loll ^ XL m,n, \j,,v lu 1/ ^ 1 OU I 

m, n — 1) /*, i/ + 1 m.,n + l,/j,,i/ — 1 

and zero everywhere else. The vector / is constructed from the elements that provide an 
independent term for the equations, that is, i? ro|Ihft „ , while the matrix M corresponds 

0,0,0,0 

to all other elements. 
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